---
title: "Price robertson covariance and sex-specific traits"
author: "Manas Geeta Arun"
date: "2024-01-17"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "C:/Users/msamant/Documents/GitHub/Sex-specific-Price-Robertson-Covariances")
```


## Load packages

```{r message=FALSE, warning=FALSE}
library(standardize)
library(MCMCglmm)
library(ggplot2)
library(cowplot)
library(lmerTest)
library(latex2exp)

diagnostics = TRUE # Should autocorrelation and time series plots for MCMCglmm models be printed?

# Load data

d = read.csv("Comprehensive_traits.csv", header = T)

```

We investigate sex and sex ratio specific selection operating on a large number of traits measured in a laboratory population of Drosophila melanogaster. First, we ask if experimentally modulating the adult sex ratio affects male and female additive genetic variances for relative fitness (Vw) as well as the intersexual additive genetic covariance for relative fitness (COVmfw). If selection gets stronger with increasing male bias in sex ratios, we should observe a corresponding increase in Vw as we move from female biased to equal to male biased sex ratio. Additionally comparing COVmfw to Vw in males and females also allows us to compare the relative magnitude of effects of selection operating through the opposite sex relative to the magnitude of the within-sex effects of selection.

# Fitness analyses (Figure 1, Table S2)

```{r echo=TRUE, va_fitness}
set.seed(1234)

d_fit = d[d$Trait=="Fitness_Female_Female_biased"|d$Trait=="Fitness_Female_Equal"|d$Trait=="Fitness_Female_Male_biased"|d$Trait=="Fitness_Male_Female_biased"|d$Trait=="Fitness_Male_Equal"|d$Trait=="Fitness_Male_Male_biased",]

# Add columns for sex and sex ratio

d_fit$Sex = ifelse(grepl("Fitness_Female", d_fit$Trait), "Female", "Male")
d_fit$Sex.Ratio = ifelse(grepl("Equal", d_fit$Trait), "Equal", ifelse(grepl("Female_biased", d_fit$Trait), "Female_biased", "Male_biased"))

d_fit$Treatment = paste(d_fit$Sex, d_fit$Sex.Ratio)
d_fit$RelFitness = d_fit$Measurement/ave(d_fit$Measurement, d_fit$Treatment)

d_fit$Sex.Ratio = as.factor(d_fit$Sex.Ratio)
d_fit$Family = as.character(d_fit$Family)
d_fit$Family = as.factor(d_fit$Family)
d_fit$Replicate = as.factor(d_fit$Replicate)


# Fit the model

prior_nd_1 = list(R = list(V = diag(6), nu = 0.002), G = list(G1 = list(V = diag(2), nu = 2, alpha.mu = c(0,0), alpha.V = diag(2)*1000), G2 = list(V = diag(2), nu = 2, alpha.mu = c(0,0), alpha.V = diag(2)*1000), G3 = list(V = diag(2), nu = 2, alpha.mu = c(0,0), alpha.V = diag(2)*1000), G4 = list(V = 1, nu = 0.002)))

fit_nd_1 <- MCMCglmm(RelFitness ~ Sex + Sex.Ratio + Sex:Sex.Ratio - 1, random = ~ us(Sex:at.level(Sex.Ratio, "Equal")):Family + us(Sex:at.level(Sex.Ratio, "Male_biased")):Family + us(Sex:at.level(Sex.Ratio, "Female_biased")):Family + Replicate:Family, rcov = ~idh(Sex:Sex.Ratio):units,family = "gaussian", nitt = 100000, burnin = 25000, thin=50, data = d_fit, prior = prior_nd_1, pr=T, verbose = F)



### Additive genetic variance

# Female

va_fem_e = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexFemale:at.level(Sex.Ratio, "Equal").Family'])
va_fem_m = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family'])
va_fem_f = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family'])

va_fem_e_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexFemale:at.level(Sex.Ratio, "Equal").Family'])
va_fem_m_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family'])
va_fem_f_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family'])



# Male

va_male_e = mean(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
va_male_m = mean(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
va_male_f = mean(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])

va_male_e_ci = HPDinterval(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
va_male_m_ci = HPDinterval(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
va_male_f_ci = HPDinterval(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])


### Additive genetic covariance between males and females

cov_mf_e = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
cov_mf_m = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
cov_mf_f = mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])


cov_mf_e_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
cov_mf_m_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
cov_mf_f_ci = HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])


### Variance associated with day:family

day_family = mean(fit_nd_1$VCV[,"Replicate:Family"])
day_family_ci = HPDinterval(fit_nd_1$VCV[,"Replicate:Family"])

### Residual variances

# Female

vr_fem_e = mean(fit_nd_1$VCV[,'SexFemale:Sex.RatioEqual.units'])
vr_fem_f = mean(fit_nd_1$VCV[,'SexFemale:Sex.RatioFemale_biased.units'])
vr_fem_m = mean(fit_nd_1$VCV[,'SexFemale:Sex.RatioMale_biased.units'])

vr_fem_e_ci = HPDinterval(fit_nd_1$VCV[,'SexFemale:Sex.RatioEqual.units'])
vr_fem_f_ci = HPDinterval(fit_nd_1$VCV[,'SexFemale:Sex.RatioFemale_biased.units'])
vr_fem_m_ci = HPDinterval(fit_nd_1$VCV[,'SexFemale:Sex.RatioMale_biased.units'])



# Male

vr_male_e = mean(fit_nd_1$VCV[,'SexMale:Sex.RatioEqual.units'])
vr_male_f = mean(fit_nd_1$VCV[,'SexMale:Sex.RatioFemale_biased.units'])
vr_male_m = mean(fit_nd_1$VCV[,'SexMale:Sex.RatioMale_biased.units'])

vr_male_e_ci = HPDinterval(fit_nd_1$VCV[,'SexMale:Sex.RatioEqual.units'])
vr_male_f_ci = HPDinterval(fit_nd_1$VCV[,'SexMale:Sex.RatioFemale_biased.units'])
vr_male_m_ci = HPDinterval(fit_nd_1$VCV[,'SexMale:Sex.RatioMale_biased.units'])




### Test for differences between male biased and female biased sex ratios for intersexual additive genetic covariance

mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])

HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])


### Test for differences between male biased and female biased sex ratios for Vw female

mean(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family'])


HPDinterval(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family'])

### Test for differences between male biased and female biased sex ratios for Vw male

mean(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])

HPDinterval(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'] - 2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])



############################
##### Plot (Figure 1) ######
############################

# Create a new data frame

data_va = data.frame("Sex" = c("Female Vw", "Female Vw", "Female Vw", "Male Vw", "Male Vw", "Male Vw", "COVmf,w", "COVmf,w", "COVmf,w"), 
                     "Sex Ratio" = c("Female Biased", "Equal", "Male Biased", "Female Biased", "Equal", "Male Biased", "Female Biased", "Equal", "Male Biased"), 
                     "Va" = round(c(va_fem_f, va_fem_e, va_fem_m, va_male_f, va_male_e, va_male_m, cov_mf_f, cov_mf_e, cov_mf_m), digits=4), 
                     "Va_L" = round(c(va_fem_f_ci[1], va_fem_e_ci[1], va_fem_m_ci[1], va_male_f_ci[1], va_male_e_ci[1], va_male_m_ci[1], cov_mf_f_ci[1], cov_mf_e_ci[1], cov_mf_m_ci[1]), digits=4), 
                     "Va_U" = round(c(va_fem_f_ci[2], va_fem_e_ci[2], va_fem_m_ci[2], va_male_f_ci[2], va_male_e_ci[2], va_male_m_ci[2], cov_mf_f_ci[2], cov_mf_e_ci[2], cov_mf_m_ci[2]), digits=4),
                     "Vr" = round(c(vr_fem_f, vr_fem_e, vr_fem_m, vr_male_f, vr_male_e, vr_male_m, NA, NA, NA), digits=4), 
                     "Vr_L" = round(c(vr_fem_f_ci[1], vr_fem_e_ci[1], vr_fem_m_ci[1], vr_male_f_ci[1], vr_male_e_ci[1], vr_male_m_ci[1], NA, NA, NA), digits=4), 
                     "Vr_U" = round(c(vr_fem_f_ci[2], vr_fem_e_ci[2], vr_fem_m_ci[2], vr_male_f_ci[2], vr_male_e_ci[2], vr_male_m_ci[2], NA, NA, NA), digits=4))
data_va$Sex.Ratio = factor(data_va$Sex.Ratio, levels = c("Female Biased", "Equal", "Male Biased"))
data_va$Sex = factor(data_va$Sex, levels = c("Female Vw", "Male Vw", "COVmf,w"))

plot_va_fitness = ggplot(data_va, aes(y = Va, x = Sex, fill = Sex.Ratio))
plot_va_fitness = plot_va_fitness + geom_bar(stat="identity", color="black", position=position_dodge()) + theme_bw() + 
  geom_errorbar(aes(ymin=Va_L, ymax=Va_U), width=.2, position=position_dodge(.9)) + 
  theme(text = element_text(size = 20)) + 
  labs(x = "", y = "Additive genetic (co)variance for relative fitness", fill = "Sex Ratio") + 
  scale_fill_manual(values = c("#E5E5E5", "#999999", "black")) +
  scale_x_discrete(labels=c("Female Vw" = TeX(r"(Female $V_{a,W}$)"), "Male Vw" = TeX(r"(Male $V_{a,W}$)"),"COVmf,w" = TeX(r"($COV_{mf,a,W}$)")))

plot_va_fitness
ggsave(plot_va_fitness,
       height = 7,
       width = 8.5,
       file = "Output/Figure_1.jpg")

###########################
#### Output (Table S2) ####
###########################

# Day:Family effects

day_family
day_family_ci

#print(data_va)


#########################
### Model diagnostics ###
#########################

if(diagnostics){
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexFemale:at.level(Sex.Ratio, "Equal").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexFemale:at.level(Sex.Ratio, "Equal").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexFemale:at.level(Sex.Ratio, "Male_biased").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexFemale:at.level(Sex.Ratio, "Female_biased").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexMale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Equal"):SexMale:at.level(Sex.Ratio, "Equal").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Male_biased"):SexMale:at.level(Sex.Ratio, "Male_biased").Family']))
  
  plot(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family'])
  print(autocorr(2*fit_nd_1$VCV[,'SexFemale:at.level(Sex.Ratio, "Female_biased"):SexMale:at.level(Sex.Ratio, "Female_biased").Family']))

  
}


```

# Testing for the significance of additive genetic variance for fitness (Figure S1)

Since estimates of variances obtained using MCMCglmm are constrained to be greater than 1 the 95% credible intervals are also constrained to not include 0. To test whether there was significant additive genetic variance for relative fitness, we used a rnadomisation based approach following Walter et al (2018, Am. Nat.). Specifically, for each randomisation, we permute the Family column of d_fit randomly and then fit the model fit_nd_1. For each randomisation, we record the posterior means for sex and sex ratio specific additive genetic and residual variances for relative fitness. This allows us to construct null distributions for these quantities.These null-distributions were constructed using the Rscript "Randomise_fitness.R" and can be found in "null_distributions/fitness_randomised_NULL.csv". We then caclulate the 95% credible intervals of these null distributions and examine if estimates of sex and sex ratio specific additive genetic variace for relative fitness fall outside.

```{r va_fitness_signific}

# Load the null distributions for sex and sex ratio specific fitness

d_null_fitness = read.csv("null_distributions/fitness_randomised_NULL.csv", header=T)

# Store 95% credible intervals for the null distributions of additive genetic (co)variances

va_fem_f_ran_ci = c(quantile(d_null_fitness$va_fem_f_ran, 0.025), quantile(d_null_fitness$va_fem_f_ran, 0.975))
va_fem_e_ran_ci = c(quantile(d_null_fitness$va_fem_e_ran, 0.025), quantile(d_null_fitness$va_fem_e_ran, 0.975))
va_fem_m_ran_ci = c(quantile(d_null_fitness$va_fem_m_ran, 0.025), quantile(d_null_fitness$va_fem_m_ran, 0.975))

va_male_f_ran_ci = c(quantile(d_null_fitness$va_male_f_ran, 0.025), quantile(d_null_fitness$va_male_f_ran, 0.975))
va_male_e_ran_ci = c(quantile(d_null_fitness$va_male_e_ran, 0.025), quantile(d_null_fitness$va_male_e_ran, 0.975))
va_male_m_ran_ci = c(quantile(d_null_fitness$va_male_m_ran, 0.025), quantile(d_null_fitness$va_male_m_ran, 0.975))

cov_mf_f_ran_ci = c(quantile(d_null_fitness$cov_mf_f_ran, 0.025), quantile(d_null_fitness$cov_mf_f_ran, 0.975))
cov_mf_e_ran_ci = c(quantile(d_null_fitness$cov_mf_e_ran, 0.025), quantile(d_null_fitness$cov_mf_e_ran, 0.975))
cov_mf_m_ran_ci = c(quantile(d_null_fitness$cov_mf_e_ran, 0.025), quantile(d_null_fitness$cov_mf_m_ran, 0.975))

# Add columns corresponding to the 95% CIs of the null distribution to data_va



data_va$Va_null_L = round(c(va_fem_f_ran_ci[1], va_fem_e_ran_ci[1], va_fem_m_ran_ci[1], va_male_f_ran_ci[1], va_male_e_ran_ci[1], va_male_m_ran_ci[1], cov_mf_f_ran_ci[1], cov_mf_e_ran_ci[1], cov_mf_m_ran_ci[1]), digits=4)
data_va$Va_null_U = round(c(va_fem_f_ran_ci[2], va_fem_e_ran_ci[2], va_fem_m_ran_ci[2], va_male_f_ran_ci[2], va_male_e_ran_ci[2], va_male_m_ran_ci[2], cov_mf_f_ran_ci[2], cov_mf_e_ran_ci[2], cov_mf_m_ran_ci[2]), digits=4)


#################
### Figure S1 ###
#################


plot_fitness_null = ggplot(data_va, aes(y = Va, x = Sex, color = Sex.Ratio))
plot_fitness_null = plot_fitness_null + 
                    geom_point(stat="identity", position=position_dodge(0.9), size = 3) + theme_bw() + 
                    geom_errorbar(aes(ymin=Va_null_L, ymax=Va_null_U), width=.3, position=position_dodge(.9), size =1) + 
                    theme(text = element_text(size = 20)) + 
                    labs(x = "", y = "Additive genetic (co)variance for relative fitness", color = "Sex Ratio") + scale_color_manual(values = c("#CCCCCC", "#999999", "black")) +
  scale_x_discrete(labels=c("Female Vw" = TeX(r"(Female $V_{a,W}$)"), "Male Vw" = TeX(r"(Male $V_{a,W}$)"),"COVmf,w" = TeX(r"($COV_{mf,a,W}$)")))

plot_fitness_null
ggsave(plot_fitness_null,
       height = 7,
       width = 8.5,
       file = "Output/Figure_S1.jpg")

# For table 1, combine Va, Va_L, and Va_U in one column with brackets

data_va$Va_for_table1 = paste(data_va$Va, " [", data_va$Va_L, ",", data_va$Va_U, "]", sep="")

data_va$Vr_for_table1 = paste(data_va$Vr, " [", data_va$Vr_L, ",", data_va$Vr_U, "]", sep="")


print(data_va)

write.csv(data_va, file = "Output/fitness_vA_output.csv", row.names = F)

```
# Examining repeatability of fitness measurements

We had performed two replicate assays for sex and sex ratio specific fitness. These replicate assays were performed using a separate set of flies generated in a different generation. Here we measure the among-replicate correlation for line means for sex and sex ratio specific fitness.

```{r repeatability_fitness_measures}

d_fit1 = aggregate(d_fit, Measurement ~ Sex*Sex.Ratio*Replicate*Family, FUN="mean")

fit_corr = function(sex, sr){
  
  d_fit2 = d_fit1[d_fit1$Sex==sex,]
  d_fit3 = d_fit2[d_fit2$Sex.Ratio==sr,]
  
  days = unique(d_fit3$Replicate)

  
  rep1 = rep(NA, 43)
  rep2 = rep(NA, 43)
  
  for (fam in 1:43){
    
    if(fam %in% d_fit3[d_fit3$Replicate==days[1],]$Family){
      rep1[fam] = d_fit3[d_fit3$Replicate==days[1]&d_fit3$Family==fam,]$Measurement
    }else{rep1[fam] = NA}
    
    if(fam %in% d_fit3[d_fit3$Replicate==days[2],]$Family){
      rep2[fam] = d_fit3[d_fit3$Replicate==days[2]&d_fit3$Family==fam,]$Measurement
    }else{rep1[fam] = NA}
  }
  
  

  plot(rep1/mean(rep1, na.rm=TRUE), rep2/mean(rep2, na.rm=TRUE), title(paste(sex, "-", sr)))
  abline(0,1)
  return(cor(rep1, rep2, use = "complete.obs"))
  
}

fit_corr("Male", "Male_biased")
fit_corr("Male", "Female_biased")
fit_corr("Male", "Equal")

fit_corr("Female", "Male_biased")
fit_corr("Female", "Female_biased")
fit_corr("Female", "Equal")

```


# Testing for additive genetic variance for traits (Figure S2, Table S5)

We used a similar approach to test if there was significant additive genetic variance for various male and female traits.We generated null distributions for additive genetic variance for traits using the Rscript "randomise_traits.R". This scripts does the following:

1. Load the comprehensive data set and select the data for fitness
2. Select data for the correct trait
3. Create a new "randomised" data set by randomly permuting the elements of the "Family" column
4. Fit model and store mean of posterior distribution
5. Repeat steps 1-4 1000 times for each trait

We then fit MCMCglmm models to our actual data and examine if actual estimates (means of posterior distributions) of Va for each trait lie outside the 95% credible intervals of the corresponding null distribution obtained by randomising data sets.

```{r va_traits}
# Load data

std = TRUE # Should traits data be standardised to have a mean 0 and stdev 1? (depends on whether using the null distribution that was generated using standardised data)

d = read.csv("Comprehensive_traits.csv", header = T)
d$Trait = factor(d$Trait)
# Make a list of all the traits  
list_traits = levels(d$Trait)

# Remove the names of fitness data.
list_traits = list_traits[grep("Fitness", list_traits, invert = T)]

# Load null distributions

d_null_trait = read.csv("null_distributions/trait_standardised_randomised_NULL.csv", header=T)

# If this dataset doesn't have a heritability column add it
# Ideally we should have a posterior mean estimate of heritability computed separately for each replicate randomisation
# If that's not done calculate the heritability from posterior means of va and vr

if(!("h2_trait"%in%names(d_null_trait))){
  d_null_trait$h2_trait = 2*d_null_trait$va_trait/(d_null_trait$va_trait + d_null_trait$vr_trait)
}

# Create empty matrix to store data

trait_variation_data = matrix(NA, nrow = length(list_traits), ncol = 21)

# Loop through the traits

count = 1

for (trait in list_traits){
  
  print(trait)
  
  # Select data
  d_trait = d[d$Trait==trait,]
  
  # Standardise the data if required
  
  if(std){d_trait$Measurement = scale(d_trait$Measurement)}
  
  # Fit model
  
  prior_trait = list(R = list(V = 1, nu = 0.002), G = list(G1 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 1000), G2 = list(V = 1, nu = 0.002)))
  
  fit_trait = MCMCglmm(Measurement ~ 1, random = ~ Family + Replicate:Family, family = "gaussian", nitt = 500000, burnin = 50000, thin=300, data = d_trait, prior = prior_trait, pr=T, verbose = F)
  
  ### Store data
  
  # Actual estimates and CIs
  
  va_trait_true = mean(2*fit_trait$VCV[, "Family"])
  va_trait_true_l = HPDinterval(2*fit_trait$VCV[, "Family"])[1]
  va_trait_true_u = HPDinterval(2* fit_trait$VCV[, "Family"])[2]
  
  vr_trait_true = mean(fit_trait$VCV[, "units"])
  vr_trait_true_l = HPDinterval(fit_trait$VCV[, "units"])[1]
  vr_trait_true_u = HPDinterval(fit_trait$VCV[, "units"])[2]
  
  h2_trait_true = mean(2*fit_trait$VCV[, "Family"]/(fit_trait$VCV[, "Family"] + fit_trait$VCV[, "units"]))
  h2_trait_true_l = HPDinterval(2*fit_trait$VCV[, "Family"]/(fit_trait$VCV[, "Family"] + fit_trait$VCV[, "units"]))[1]
  h2_trait_true_u = HPDinterval(2*fit_trait$VCV[, "Family"]/(fit_trait$VCV[, "Family"] + fit_trait$VCV[, "units"]))[2]
  
  day_family_true = mean(fit_trait$VCV[, "Replicate:Family"])
  day_family_true_l = HPDinterval(fit_trait$VCV[, "Replicate:Family"])[1]
  day_family_true_u = HPDinterval(fit_trait$VCV[, "Replicate:Family"])[2]
  
  # CIs from null distribution
  
  va_trait_null_l = quantile(d_null_trait[d_null_trait$Trait==trait,]$va_trait, 0.025)
  va_trait_null_u = quantile(d_null_trait[d_null_trait$Trait==trait,]$va_trait, 0.975)
  
  vr_trait_null_l = quantile(d_null_trait[d_null_trait$Trait==trait,]$vr_trait, 0.025)
  vr_trait_null_u = quantile(d_null_trait[d_null_trait$Trait==trait,]$vr_trait, 0.975)
  
  h2_trait_null_l = quantile(d_null_trait[d_null_trait$Trait==trait,]$h2_trait, 0.025)
  h2_trait_null_u = quantile(d_null_trait[d_null_trait$Trait==trait,]$h2_trait, 0.975)
  
  day_family_null_l = quantile(d_null_trait[d_null_trait$Trait==trait,]$day_family_ran, 0.025)
  day_family_null_u = quantile(d_null_trait[d_null_trait$Trait==trait,]$day_family_ran, 0.975)
  
  # Model diagnostics
  
  if(diagnostics){
    plot(2*fit_trait$VCV[, "Family"])
    print(autocorr(2*fit_trait$VCV[, "Family"]))
  }
  
  # Add data of the current trait to the data frame
  
  trait_variation_data[count,] = c(trait,
                                   round(va_trait_true, digits=4), round(va_trait_true_l, digits=4), round(va_trait_true_u, digits=4),
                                   round(vr_trait_true, digits=4), round(vr_trait_true_l, digits=4), round(vr_trait_true_u, digits=4),
                                   round(h2_trait_true, digits=4), round(h2_trait_true_l, digits=4), round(h2_trait_true_u, digits=4),
                                   round(day_family_true, digits=4), round(day_family_true_l, digits=4), round(day_family_true_u, digits=4),
                                   round(va_trait_null_l, digits=4), round(va_trait_null_u, digits=4),
                                   round(vr_trait_null_l, digits=4), round(vr_trait_null_u, digits=4),
                                   round(h2_trait_null_l, digits=4), round(h2_trait_null_u, digits=4),
                                   round(day_family_null_l, digits=4), round(day_family_null_u, digits=4))
  
  count = count + 1

  
}

colnames(trait_variation_data) = c("trait",
                                   "va_trait_true", "va_trait_true_l", "va_trait_true_u",
                                   "vr_trait_true", "vr_trait_true_l", "vr_trait_true_u",
                                   "h2_trait_true", "h2_trait_true_l", "h2_trait_true_u",
                                   "day_family_true", "day_family_true_l", "day_family_true_u",
                                   "va_trait_null_l", "va_trait_null_u",
                                   "vr_trait_null_l", "vr_trait_null_u",
                                   "h2_trait_null_l", "h2_trait_null_u",
                                   "day_family_null_l", "day_family_null_u")
trait_variation_data = data.frame(trait_variation_data)

print(trait_variation_data)

# Add columns that merge point estimates and CIs together (helps create Tables for the manuscript)

trait_variation_data$trait_names_for_table = c('Copulation duration in females','Copulation duration in males','Female fecundity (female biased)', 'Female fecundity (male biased)', 'Male effect on female fecundity (female biased)', 'Male effect on female fecundity (male biased)', 'Female latency to first mating', 'Male latency to first mating', 'Female mating rate (female biased)', 'Female mating rate (male biased)', 'Male mating rate (female biased)', 'Male mating rate (male biased)', 'Sperm defence ability (P1)', 'Sperm offence ability (P2)', 'Female remating latency', 'Male mating latency with mated females')

trait_variation_data$va_for_table = paste(trait_variation_data$va_trait_true, " [", trait_variation_data$va_trait_true_l, ",", trait_variation_data$va_trait_true_u, "]", sep = "")

trait_variation_data$vr_for_table = paste(trait_variation_data$vr_trait_true, " [", trait_variation_data$vr_trait_true_l, ",", trait_variation_data$vr_trait_true_u, "]", sep = "")

trait_variation_data$day_family_for_table = paste(trait_variation_data$day_family_true, " [", trait_variation_data$day_family_true_l, ",", trait_variation_data$day_family_true_u, "]", sep = "")

write.csv(trait_variation_data, file = "Output/trait_vA_output.csv", row.names = F)

#################
### Figure S2 ###
#################

# Make plot of va

p_trait_va = ggplot(trait_variation_data, aes(y = as.numeric(va_trait_true), x = trait))
p_trait_va = p_trait_va + geom_point(size = 2) + theme_bw() + 
  geom_errorbar(aes(ymin=as.numeric(va_trait_null_l), ymax=as.numeric(va_trait_null_u)), width=.5, size =0.5, position=position_dodge(1)) + 
  theme(text = element_text(size = 15)) + 
  labs(x = "Trait", y = "Additive genetic variance on standardised scale") + 
  scale_x_discrete(labels = c('Copulation duration in females','Copulation duration in males','Female fecundity (female biased)', 'Female fecundity (male biased)', 'Male effect on female fecundity (female biased)', 'Male effect on female fecundity (male biased)', 'Female latency to first mating', 'Male latency to first mating', 'Female mating rate (female biased)', 'Female mating rate (male biased)', 'Male mating rate (female biased)', 'Male mating rate (male biased)', 'Sperm defence ability (P1)', 'Sperm offence ability (P2)', 'Female remating latency', 'Male mating latency with mated females')) +
  coord_flip() 

p_trait_va
ggsave(p_trait_va,
       height = 5,
       width = 8.5, 
       file = "Output/Figure_S2.jpg")

  

```

# Traits (Figures 2-4, Tables S3-S4)

Next, for each of the traits we measured (irrespective of whether the trait was measured in males or females), we investigated the additive genetic covariance between the trait and either male fitness (COVm) (at male biased, equal and female biased sex ratios) or female fitness (at male biased, equal and female biased sex ratios) (COVf). Sex-specific Price-Robertson Identity predicts that the expected change in the trait mean due to one generation of selection is equal to the mean of COVm and COVf.

Plots are generated within "prc_analysis.html", while numerical output is stored within "PRC_output.csv" (for "non-sex ratio traits") and "PRC_output_sep_sr.csv" (for "sex ratio traits")


## Function to calculate Price-Robertson Covariance

This function first identifies whether the trait is a "sex-ratio trait" or a "non-sex ratio trait" (see Statistical Analysis)
- Computes within-sex and across-sex PRCs
- Prepares plots
- Stores outputs (estimates of PRCs, estimates of heritabilities, etc.) within "Output/PRC_output.csv" (non-sex ratio traits) and "Output/PRC_output_sep_sr.csv" (sex-ratio traits)


```{r prc_function}
prc = function(trait, # Name of the trait
               diagnostics=FALSE, # print autocorrelation and time series plots? 
               d=d ){   # Comprehensive trait data)
  
  # Convert trait into a factor
  
  d$Trait = factor(d$Trait)
  
  # Divide fitness data by means to obtain relative fitnesses
  
  d[d$Trait=="Fitness_Female_Equal",]$Measurement = (d[d$Trait=="Fitness_Female_Equal",]$Measurement)/(mean(d[d$Trait=="Fitness_Female_Equal",]$Measurement))
  d[d$Trait=="Fitness_Female_Female_biased",]$Measurement = (d[d$Trait=="Fitness_Female_Female_biased",]$Measurement)/(mean(d[d$Trait=="Fitness_Female_Female_biased",]$Measurement))
  d[d$Trait=="Fitness_Female_Male_biased",]$Measurement = (d[d$Trait=="Fitness_Female_Male_biased",]$Measurement)/(mean(d[d$Trait=="Fitness_Female_Male_biased",]$Measurement))
  d[d$Trait=="Fitness_Male_Male_biased",]$Measurement = (d[d$Trait=="Fitness_Male_Male_biased",]$Measurement)/(mean(d[d$Trait=="Fitness_Male_Male_biased",]$Measurement))
  d[d$Trait=="Fitness_Male_Female_biased",]$Measurement = (d[d$Trait=="Fitness_Male_Female_biased",]$Measurement)/(mean(d[d$Trait=="Fitness_Male_Female_biased",]$Measurement))
  d[d$Trait=="Fitness_Male_Equal",]$Measurement = (d[d$Trait=="Fitness_Male_Equal",]$Measurement)/(mean(d[d$Trait=="Fitness_Male_Equal",]$Measurement))
  
  #### Identify if the trait is measured in a sex-ratio environment. If yes only analyse with two sex ratios (sr = 2), else sr = 3
  
  # List of traits measured at sex ratio environments (also exists prc_analysis.Rmd. If editing here, also edit there!!!)
  
  sr_list = c("Female_fecundity","Mate_fecundity" , "MR_Male", "MR_Female", "LA_Female", "LA_Male", "LA_P1_Female", "LA_P1_Male", "LA_P2_Female", "LA_P2_Male", "LA_P3_Female", "LA_P3_Male")
  
  sr = ifelse(trait %in% sr_list, 2, 3)
  
  if(sr==3){
    
    # Standardise the trait data
    
    d[d$Trait==trait,]$Measurement = scale(d[d$Trait==trait,]$Measurement)
    
    
    
    ### trait and male fitness ###
    
    # Create a separate data frame containing the data for trait and male fitness
    
    d_maleW = d[d$Trait==trait|d$Trait=="Fitness_Male_Equal"|d$Trait=="Fitness_Male_Male_biased"|d$Trait=="Fitness_Male_Female_biased",]
    
    set.seed(15081947)
    
    prior_maleW = list(R = list(V = diag(4), nu = 0.002), G = list(G1 = list(V = diag(4), nu = 2, alpha.mu = c(0,0,0,0), alpha.V = diag(4)*1000), G2 = list(V = 1, nu = 0.002)))
    
    fit_maleW = MCMCglmm(Measurement ~ Trait -1, random = ~us(Trait):Family + Replicate:Trait, rcov = ~idh(Trait):units, family = "gaussian", nitt = 100000, burnin = 25000, thin=50, data = d_maleW, prior = prior_maleW, pr=T, verbose = F)
    
    # Female biased sex ratio
    cov_trait_maleW_fb = mean(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")])
    cov_trait_maleW_fb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")])
    
    # Equal sex ratio
    cov_trait_maleW_e = mean(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Equal.Family", sep = "")])
    cov_trait_maleW_e_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Equal.Family", sep = "")])
    
    # Male biased sex ratio
    cov_trait_maleW_mb = mean(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")])
    cov_trait_maleW_mb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")])
    
    # Difference between male biased and female biased sex ratios
    diff_trait_maleW = mean(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")] - 2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")])
    diff_trait_maleW_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")] - 2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")])
    
    ## Trait heritability
    
    h2_trait_m = mean(2*fit_maleW$VCV[,paste("Trait", trait, ":Trait", trait, ".Family", sep = "")]) # The "_m" indicates that this is an estimate from the model which also had female fitness data
    h2_trait_m_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, ":Trait", trait, ".Family", sep = "")])
    
    ## Replicate effects
    
    rep_m = mean(fit_maleW$VCV[, "Replicate:Trait"])
    rep_m_ci = HPDinterval(fit_maleW$VCV[, "Replicate:Trait"])
    
    ### Diagnostics for model estimating RCs with male fitness
    
    if(diagnostics){
      
      plot(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")])
      print(autocorr(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Female_biased.Family", sep = "")]))
      
      plot(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Equal.Family", sep = "")])
      print(autocorr(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Equal.Family", sep = "")]))
      
      plot(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")])
      print(autocorr(2*fit_maleW$VCV[,paste("Trait", trait, ":TraitFitness_Male_Male_biased.Family", sep = "")]))
      
    }
    
    
    ### Trait and female fitness ###
    
    # Create a separate data frame containing the data for ML and female fitness
    
    d_femaleW = d[d$Trait==trait|d$Trait=="Fitness_Female_Equal"|d$Trait=="Fitness_Female_Male_biased"|d$Trait=="Fitness_Female_Female_biased",]
    
    set.seed(15081947)
    
    prior_femaleW = list(R = list(V = diag(4), nu = 0.002), G = list(G1 = list(V = diag(4), nu = 2, alpha.mu = c(0,0,0,0), alpha.V = diag(4)*1000), G2 = list(V = 1, nu = 0.002)))
    
    fit_femaleW = MCMCglmm(Measurement ~ Trait -1, random = ~us(Trait):Family + Replicate:Trait, rcov = ~idh(Trait):units, family = "gaussian", nitt = 100000, burnin = 25000, thin=50, data = d_femaleW, prior = prior_femaleW, pr=T, verbose = F)
    
    # Female biased sex ratio
    cov_trait_femaleW_fb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")])
    cov_trait_femaleW_fb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")])
    
    # Equal sex ratio
    cov_trait_femaleW_e = mean(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Equal.Family", sep = "")])
    cov_trait_femaleW_e_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Equal.Family", sep = "")])
    
    # Male biased sex ratio
    cov_trait_femaleW_mb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")])
    cov_trait_femaleW_mb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")])
    
    # Difference between female biased and fefemale biased sex ratios
    diff_trait_femaleW = mean(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")] - 2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")])
    diff_trait_femaleW_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")] - 2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")])
    
    ## Trait heritability
    
    h2_trait_f = mean(2*fit_femaleW$VCV[,paste("Trait", trait, ":Trait", trait, ".Family", sep = "")]) # The "_f" indicates that this is an estimate from the model which also had female fitness data
    h2_trait_f_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, ":Trait", trait, ".Family", sep = "")])
    
    ## Replicate effects
    
    rep_f = mean(fit_femaleW$VCV[, "Replicate:Trait"])
    rep_f_ci = HPDinterval(fit_femaleW$VCV[, "Replicate:Trait"])
    
    ### Diagnostics for model estimating RCs with female fitness
    
    if(diagnostics){
      
      plot(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")])
      print(autocorr(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Female_biased.Family", sep = "")]))
      
      plot(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Equal.Family", sep = "")])
      print(autocorr(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Equal.Family", sep = "")]))
      
      plot(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")])
      print(autocorr(2*fit_femaleW$VCV[,paste("Trait", trait, ":TraitFitness_Female_Male_biased.Family", sep = "")]))
      
    }
    
    
    ################# PLOT #####################
    
    # Create a new data frame
    
    data_plot_trait = data.frame("Sex" = c("Female", "Female", "Female", "Male", "Male", "Male"), "Sex Ratio" = c("Female Biased", "Equal", "Male Biased", "Female Biased", "Equal", "Male Biased"), "Cov" = c(cov_trait_femaleW_fb, cov_trait_femaleW_e, cov_trait_femaleW_mb, cov_trait_maleW_fb, cov_trait_maleW_e, cov_trait_maleW_mb), "Cov_L" = c(cov_trait_femaleW_fb_ci[1], cov_trait_femaleW_e_ci[1], cov_trait_femaleW_mb_ci[1], cov_trait_maleW_fb_ci[1], cov_trait_maleW_e_ci[1], cov_trait_maleW_mb_ci[1]), "Cov_U" = c(cov_trait_femaleW_fb_ci[2], cov_trait_femaleW_e_ci[2], cov_trait_femaleW_mb_ci[2], cov_trait_maleW_fb_ci[2], cov_trait_maleW_e_ci[2], cov_trait_maleW_mb_ci[2]))
    data_plot_trait$Sex.Ratio = factor(data_plot_trait$Sex.Ratio, levels = c("Female Biased", "Equal", "Male Biased"))
    
    plot_trait = ggplot(data_plot_trait, aes(y = Cov, x = Sex, fill = Sex.Ratio)) + 
      geom_bar(stat="identity", color="black", position=position_dodge()) + theme_bw() + 
      geom_errorbar(aes(ymin=Cov_L, ymax=Cov_U), width=.2, position=position_dodge(.9)) + 
      theme(text = element_text(size = 15)) + 
      labs(y = "Robertson covariance", fill = "Sex Ratio", x = "Selected sex") + 
      scale_fill_manual(values = c("#E5E5E5", "#999999", "black")) + 
      ggtitle(trait) + 
      theme(plot.title = element_text(hjust = 0.5, size = 15))
    

    return(list(plot_trait = plot_trait,
                cov_trait_femaleW_fb = round(cov_trait_femaleW_fb, digits=4),
                cov_trait_femaleW_fb_ci = round(cov_trait_femaleW_fb_ci, digits=4),
                cov_trait_femaleW_e = round(cov_trait_femaleW_e, digits=4),
                cov_trait_femaleW_e_ci = round(cov_trait_femaleW_e_ci, digits=4),
                cov_trait_femaleW_mb = round(cov_trait_femaleW_mb, digits=4),
                cov_trait_femaleW_mb_ci = round(cov_trait_femaleW_mb_ci, digits=4),
                diff_trait_femaleW = round(diff_trait_femaleW, digits=4),
                diff_trait_femaleW_ci = round(diff_trait_femaleW_ci, digits=4),
                cov_trait_maleW_fb = round(cov_trait_maleW_fb, digits=4),
                cov_trait_maleW_fb_ci = round(cov_trait_maleW_fb_ci, digits=4),
                cov_trait_maleW_e = round(cov_trait_maleW_e, digits=4),
                cov_trait_maleW_e_ci = round(cov_trait_maleW_e_ci, digits=4),
                cov_trait_maleW_mb = round(cov_trait_maleW_mb, digits=4),
                cov_trait_maleW_mb_ci = round(cov_trait_maleW_mb_ci, digits=4),
                diff_trait_maleW = round(diff_trait_maleW, digits=4),
                diff_trait_maleW_ci = round(diff_trait_maleW_ci, digits=4),
                h2_trait_m = round(h2_trait_m, digits=4),
                h2_trait_m_ci = round(h2_trait_m_ci, digits=4),
                h2_trait_f = round(h2_trait_f, digits=4),
                h2_trait_f_ci = round(h2_trait_f_ci, digits=4),
                var_trait = round(var(d[d$Trait==trait,]$Measurement), digits=4), # Just to check that this is 1
                rep_f = round(rep_f, digits=4), # Relicate variance in model with female fitness
                rep_f_ci = round(rep_f_ci, digits=4),
                rep_m = round(rep_m, digits=4), # Replicate variance in moel with male fitness
                rep_m_ci = round(rep_m_ci, digits=4))) 
    
  }
  
  if(sr==2){
    
    # Standardise the traits
    
    d[d$Trait==paste(trait, "_Female_biased", sep = ""),]$Measurement = scale(d[d$Trait==paste(trait, "_Female_biased", sep = ""),]$Measurement)
    d[d$Trait==paste(trait, "_Male_biased", sep = ""),]$Measurement = scale(d[d$Trait==paste(trait, "_Male_biased", sep = ""),]$Measurement)
    
    
    ### Trait and male fitness ###
    
    # Create a separate data frame containing the data for male mating rate and male fitness
    
    d_maleW = d[d$Trait==paste(trait, "_Female_biased", sep = "")|d$Trait==paste(trait, "_Male_biased", sep = "")|d$Trait=="Fitness_Male_Male_biased"|d$Trait=="Fitness_Male_Female_biased",]
    
    set.seed(15081947)
    
    prior_maleW = list(R = list(V = diag(4), nu = 0.002), G = list(G1 = list(V = diag(4), nu = 2, alpha.mu = c(0,0,0,0), alpha.V = diag(4)*1000), G2 = list(V = 1, nu = 0.002)))
    
    fit_maleW = MCMCglmm(Measurement ~ Trait -1, random = ~us(Trait):Family + Replicate:Trait, rcov = ~idh(Trait):units, family = "gaussian", nitt = 100000, burnin = 25000, thin=50, data = d_maleW, prior = prior_maleW, pr=T, verbose = F)
    
    # Female biased sex ratio
    cov_trait_maleW_fb = mean(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")])
    cov_trait_maleW_fb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")])
    
    # Male biased sex ratio
    cov_trait_maleW_mb = mean(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")])
    cov_trait_maleW_mb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")])
    
    # Difference between male biased and female biased sex ratios
    diff_trait_maleW = mean(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")] - 2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")])
    diff_trait_maleW_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")] - 2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")])
    
    ### Trait heritabilities (The "_m" indicates that the estimates are from the model that also has the male fitness data)
    
    h2_trait_m_fb = mean(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":Trait", trait, "_Female_biased", ".Family", sep = "")])
    h2_trait_m_fb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":Trait", trait, "_Female_biased", ".Family", sep = "")])
    
    h2_trait_m_mb = mean(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":Trait", trait, "_Male_biased", ".Family", sep = "")])
    h2_trait_m_mb_ci = HPDinterval(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":Trait", trait, "_Male_biased", ".Family", sep = "")])
    
    ## Replicate effects
    
    rep_m = mean(fit_maleW$VCV[, "Replicate:Trait"])
    rep_m_ci = HPDinterval(fit_maleW$VCV[, "Replicate:Trait"])
    
    ### Diagnostics for model estimating RCs with male fitness
    
    if(diagnostics){
      
      plot(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")])
      print(autocorr(2*fit_maleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Male_Female_biased.Family", sep = "")]))
      
      plot(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")])
      print(autocorr(2*fit_maleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Male_Male_biased.Family", sep = "")]))
      
    }
    
    ### Trait and female fitness ###
    
    d_femaleW = d[d$Trait==paste(trait, "_Female_biased", sep = "")|d$Trait==paste(trait, "_Male_biased", sep = "")|d$Trait=="Fitness_Female_Male_biased"|d$Trait=="Fitness_Female_Female_biased",]
    
    set.seed(15081947)
    
    prior_femaleW = list(R = list(V = diag(4), nu = 0.002), G = list(G1 = list(V = diag(4), nu = 2, alpha.mu = c(0,0,0,0), alpha.V = diag(4)*1000), G2 = list(V = 1, nu = 0.002)))
    
    fit_femaleW = MCMCglmm(Measurement ~ Trait -1, random = ~us(Trait):Family + Replicate:Trait, rcov = ~idh(Trait):units, family = "gaussian", nitt = 100000, burnin = 25000, thin=50, data = d_femaleW, prior = prior_femaleW, pr=T, verbose = F)
    
    # Female biased sex ratio
    cov_trait_femaleW_fb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")])
    cov_trait_femaleW_fb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")])
    
    # Male biased sex ratio
    cov_trait_femaleW_mb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")])
    cov_trait_femaleW_mb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")])
    
    # Difference between female biased and male biased sex ratios
    diff_trait_femaleW = mean(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")] - 2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")])
    diff_trait_femaleW_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")] - 2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")])
    
    ### Trait heritabilities (The "_m" indicates that the estimates are from the model that also has the male fitness data)
    
    h2_trait_f_fb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":Trait", trait, "_Female_biased", ".Family", sep = "")])
    h2_trait_f_fb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":Trait", trait, "_Female_biased", ".Family", sep = "")])
    
    h2_trait_f_mb = mean(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":Trait", trait, "_Male_biased", ".Family", sep = "")])
    h2_trait_f_mb_ci = HPDinterval(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":Trait", trait, "_Male_biased", ".Family", sep = "")])
    
    ## Replicate effects
    
    rep_f = mean(fit_femaleW$VCV[, "Replicate:Trait"])
    rep_f_ci = HPDinterval(fit_femaleW$VCV[, "Replicate:Trait"])
    
    ### Diagnostics for model estimating RCs with female fitness
    
    if(diagnostics){
      
      plot(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")])
      print(autocorr(2*fit_femaleW$VCV[,paste("Trait", trait, "_Female_biased", ":TraitFitness_Female_Female_biased.Family", sep = "")]))
      
      plot(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")])
      print(autocorr(2*fit_femaleW$VCV[,paste("Trait", trait, "_Male_biased", ":TraitFitness_Female_Male_biased.Family", sep = "")]))
      
    }
    
    
    # Create a new data frame
    
    data_plot_trait = data.frame("Sex" = c("Female", "Female", "Male", "Male"), "Sex Ratio" = c("Female Biased", "Male Biased", "Female Biased", "Male Biased"), "Cov" = c(cov_trait_femaleW_fb, cov_trait_femaleW_mb, cov_trait_maleW_fb, cov_trait_maleW_mb), "Cov_L" = c(cov_trait_femaleW_fb_ci[1], cov_trait_femaleW_mb_ci[1], cov_trait_maleW_fb_ci[1],  cov_trait_maleW_mb_ci[1]), "Cov_U" = c(cov_trait_femaleW_fb_ci[2], cov_trait_femaleW_mb_ci[2], cov_trait_maleW_fb_ci[2], cov_trait_maleW_mb_ci[2]))
    data_plot_trait$Sex.Ratio = factor(data_plot_trait$Sex.Ratio, levels = c("Female Biased", "Male Biased"))
    
    plot_trait = ggplot(data_plot_trait, aes(y = Cov, x = Sex, fill = Sex.Ratio)) + 
      geom_bar(stat="identity", color="black", position=position_dodge()) + theme_bw() + 
      geom_errorbar(aes(ymin=Cov_L, ymax=Cov_U), width=.2, position=position_dodge(.9)) + 
      theme(text = element_text(size = 15)) + 
      labs(y = "Robertson covariance", fill = "Sex Ratio", x = "Selected sex") + 
      scale_fill_manual(values = c("#E5E5E5", "black")) + 
      ggtitle(trait) + 
      theme(plot.title = element_text(hjust = 0.5, size = 15)) 
    
    return(list(plot_trait = plot_trait,
                cov_trait_femaleW_fb = round(cov_trait_femaleW_fb, digits=4),
                cov_trait_femaleW_fb_ci = round(cov_trait_femaleW_fb_ci, digits=4),
                cov_trait_femaleW_e = NA,
                cov_trait_femaleW_e_ci = NA,
                cov_trait_femaleW_mb = round(cov_trait_femaleW_mb, digits=4),
                cov_trait_femaleW_mb_ci = round(cov_trait_femaleW_mb_ci, digits=4),
                diff_trait_femaleW = round(diff_trait_femaleW, digits=4),
                diff_trait_femaleW_ci = round(diff_trait_femaleW_ci, digits=4),
                cov_trait_maleW_fb = round(cov_trait_maleW_fb, digits=4),
                cov_trait_maleW_fb_ci = round(cov_trait_maleW_fb_ci, digits=4),
                cov_trait_maleW_e = NA, 
                cov_trait_maleW_e_ci = NA,
                cov_trait_maleW_mb = round(cov_trait_maleW_mb, digits=4),
                cov_trait_maleW_mb_ci = round(cov_trait_maleW_mb_ci, digits=4),
                diff_trait_maleW = round(diff_trait_maleW, digits=4),
                diff_trait_maleW_ci = round(diff_trait_maleW_ci, digits=4),
                h2_trait_m_fb = round(h2_trait_m_fb, digits=4),
                h2_trait_m_fb_ci = round(h2_trait_m_fb_ci, digits=4),
                h2_trait_m_mb = round(h2_trait_m_mb, digits=4),
                h2_trait_m_mb_ci = round(h2_trait_m_mb_ci, digits=4),
                h2_trait_f_fb = round(h2_trait_f_fb, digits=4),
                h2_trait_f_fb_ci = round(h2_trait_f_fb_ci, digits=4),
                h2_trait_f_mb = round(h2_trait_f_mb, digits=4),
                h2_trait_f_mb_ci = round(h2_trait_f_mb_ci, digits=4),
                var_trait_fb = round(var(d[d$Trait==paste(trait, "_Female_biased", sep = ""),]$Measurement), digits=4),
                var_trait_mb = round(var(d[d$Trait==paste(trait, "_Male_biased", sep = ""),]$Measurement), digits=4),
                rep_f = round(rep_f, digits=4), # Replicate variance in the model with female fitness
                rep_f_ci = round(rep_f_ci, digits=4),
                rep_m = round(rep_m, digits=4), # Replicate variance in the model with male fitness
                rep_m_ci = round(rep_m_ci, digits=4)))
    
  }
}


```

## Standardise trait data

```{r standardise_traits}

d$Trait = factor(d$Trait)

# Make a list of all the traits to automate implementation of "prc". 
list_traits = levels(d$Trait)

# Remove the names of fitness data.
list_traits = list_traits[grep("Fitness", list_traits, invert = T)]

# For traits measured at different sex ratios (eg, mate fecundity, mating rate) the sex-ratio suffix ("_Male_biased", etc) has to be removed for the application of prc

list_traits = gsub("_Female_biased", "", list_traits)
list_traits = gsub("_Equal", "", list_traits)
list_traits = gsub("_Male_biased", "", list_traits)

# This results in duplication of the names of traits measured at different sex ratios
# Selected only unique names

list_traits = unique(list_traits)

# Make a list of names of traits that should appear in plot titles

list_names = c("Copulation duration in females", "Copulation duration in males","Female fecundity", "Male effect on female fecundity", "Female latency to first mating", "Male latency to first mating", "Female mating rate", "Male mating rate", "Sperm defence ability (P1)", "Sperm offence ability (P2)", "Female remating latency", "Male mating latency with mated females")

print(cbind(list_names, list_traits))

```

## Plots (Panels within Figures 2-4)

Numerical outputs can be found in the files "prc_output_sep_sr.csv" and "prc_output.csv".

```{r prc_analysis_plots}
# Analyse the traits recursively

#### Create empty tables to store output
# Traits measured in sex ratio environments vs traits measured without sex ratio environments have slightly different output terms
# Therefore need separate tables for them

output_table = c() # Output table for traits measured outside sex ratio environments (eg P1, ML, CD, BW, DT, etc.)
output_table_sep = c() # Output table for traits measured at female biased or male biased sex ratios (eg. Mate_fec, LA, MR, etc.)

# List of traits measured at sex ratio environments (also exists inside prc_function!!!!! If editing here, also edit there!!!)

sr_list = c("Female_fecundity","Mate_fecundity" , "MR_Male", "MR_Female")

# Create an empty list to store plot objects
plot_list = list()
count = 1 # Just for indexing this list in the following loop

for (trait in list_traits){
  
  print(trait)
  prc_trait = prc(trait=trait, diagnostics=diagnostics, d=d)
  
  # Print plot
  #print(prc_trait$plot_trait + ggtitle(list_names[which(list_traits == trait)]))
  
  plot_list[[count]] = prc_trait$plot_trait + ggtitle(list_names[which(list_traits == trait)])
  count = count + 1
  
  ## Compile the output in a table
  # Traits measured in sex ratio environments vs traits measured without sex ratio environments have slightly different output terms
  # Therefore need separate tables for them
  
  if (trait %in% sr_list){
    
      current_output = c(trait, 
                       paste(prc_trait$cov_trait_femaleW_fb, " [", prc_trait$cov_trait_femaleW_fb_ci[1], ", ", prc_trait$cov_trait_femaleW_fb_ci[2], "]", sep = ""), 
                       paste(prc_trait$cov_trait_femaleW_mb, " [", prc_trait$cov_trait_femaleW_mb_ci[1], ", ", prc_trait$cov_trait_femaleW_mb_ci[2], "]", sep = ""), 
                       paste(prc_trait$diff_trait_femaleW, " [", prc_trait$diff_trait_femaleW_ci[1], ", ", prc_trait$diff_trait_femaleW_ci[2], "]", sep = ""), 
                       paste(prc_trait$cov_trait_maleW_fb, " [", prc_trait$cov_trait_maleW_fb_ci[1], ", ",  prc_trait$cov_trait_maleW_fb_ci[2], "]", sep =""),
                       paste(prc_trait$cov_trait_maleW_mb, " [", prc_trait$cov_trait_maleW_mb_ci[1], ", ",  prc_trait$cov_trait_maleW_mb_ci[2], "]", sep =""), 
                       paste(prc_trait$diff_trait_maleW, " [", prc_trait$diff_trait_maleW_ci[1], ", ", prc_trait$diff_trait_maleW_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_f_fb, " [", prc_trait$h2_trait_f_fb_ci[1], ", ", prc_trait$h2_trait_f_fb_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_f_mb, " [", prc_trait$h2_trait_f_mb_ci[1], ", ", prc_trait$h2_trait_f_mb_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_m_fb, " [", prc_trait$h2_trait_m_fb_ci[1], ", ", prc_trait$h2_trait_m_fb_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_m_mb, " [", prc_trait$h2_trait_m_mb_ci[1], ", ", prc_trait$h2_trait_m_mb_ci[2], "]", sep = ""),
                       paste(prc_trait$rep_f, " [", prc_trait$rep_f_ci[1], ", ", prc_trait$rep_f_ci[2], "]", sep = ""),
                       paste(prc_trait$rep_m, " [", prc_trait$rep_m_ci[1], ", ", prc_trait$rep_m_ci[2], "]", sep = ""))
    
    output_table_sep = rbind(output_table_sep, current_output)
    
  }else{
      
    current_output = c(trait, 
                       paste(prc_trait$cov_trait_femaleW_fb, " [", prc_trait$cov_trait_femaleW_fb_ci[1], ", ", prc_trait$cov_trait_femaleW_fb_ci[2], "]", sep = ""), 
                       paste(prc_trait$cov_trait_femaleW_e, " [", prc_trait$cov_trait_femaleW_e_ci[1], ", ", prc_trait$cov_trait_femaleW_e_ci[2], "]", sep = ""), 
                       paste(prc_trait$cov_trait_femaleW_mb, " [", prc_trait$cov_trait_femaleW_mb_ci[1], ", ", prc_trait$cov_trait_femaleW_mb_ci[2], "]", sep = ""), 
                       paste(prc_trait$diff_trait_femaleW, " [", prc_trait$diff_trait_femaleW_ci[1], ", ", prc_trait$diff_trait_femaleW_ci[2], "]", sep = ""), 
                       paste(prc_trait$cov_trait_maleW_fb, " [", prc_trait$cov_trait_maleW_fb_ci[1], ", ",  prc_trait$cov_trait_maleW_fb_ci[2], "]", sep =""),
                       paste(prc_trait$cov_trait_maleW_e, " [", prc_trait$cov_trait_maleW_e_ci[1], ", ",  prc_trait$cov_trait_maleW_e_ci[2], "]", sep =""), 
                       paste(prc_trait$cov_trait_maleW_mb, " [", prc_trait$cov_trait_maleW_mb_ci[1], ", ",  prc_trait$cov_trait_maleW_mb_ci[2], "]", sep =""), 
                       paste(prc_trait$diff_trait_maleW, " [", prc_trait$diff_trait_maleW_ci[1], ", ", prc_trait$diff_trait_maleW_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_f, " [", prc_trait$h2_trait_f_ci[1], ", ", prc_trait$h2_trait_f_ci[2], "]", sep = ""),
                       paste(prc_trait$h2_trait_m, " [", prc_trait$h2_trait_m_ci[1], ", ", prc_trait$h2_trait_m_ci[2], "]", sep = ""),
                       paste(prc_trait$rep_f, " [", prc_trait$rep_f_ci[1], ", ", prc_trait$rep_f_ci[2], "]", sep = ""),
                       paste(prc_trait$rep_m, " [", prc_trait$rep_m_ci[1], ", ", prc_trait$rep_m_ci[2], "]", sep = ""))
    
    output_table = rbind(output_table, current_output)
    
  }
  
 


}

# Format the tables and print

colnames(output_table) = c("trait", "cov_trait_femaleW_fb", "cov_trait_femaleW_e", "cov_trait_femaleW_mb", "diff_trait_femaleW", "cov_trait_maleW_fb", "cov_trait_maleW_e", "cov_trait_maleW_mb", "diff_trait_maleW", "h2_trait_f", "h2_trait_m", "rep_f", "rep_m")
rownames(output_table) = c()

colnames(output_table_sep) = c("trait", "cov_trait_femaleW_fb", "cov_trait_femaleW_mb", "diff_trait_femaleW", "cov_trait_maleW_fb", "cov_trait_maleW_mb", "diff_trait_maleW", "h2_trait_f_fb", "h2_trait_f_mb", "h2_trait_m_fb", "h2_trait_m_mb", "rep_f", "rep_m")
rownames(output_table_sep) = c()

# Output for data presented in Tables S3-S4 (PRCs), and Tables S6-S7 (heritabilities)

write.csv(output_table, "Output/PRC_output.csv", row.names = F)
write.csv(output_table_sep, "Output/PRC_output_sep_sr.csv", row.names = F)

####################
### Make Figures ###
####################

# name the list of plot objects
names(plot_list) = list_traits

# Extract legends
legend_2 = get_plot_component(plot_list$Female_fecundity + guides(color = guide_legend(nrow = 1)) + theme(legend.position = "bottom"), "guide-box-bottom") # for plots with 2 sex ratios
legend_3 = get_plot_component(plot_list$CD_female + guides(color = guide_legend(nrow = 1)) + theme(legend.position = "bottom"), "guide-box-bottom") # for plots with 3 sex ratios

# Figure 2

# Make a combined plot without legends
figure_2 = plot_grid(plot_list$Female_fecundity + theme(legend.position="none"), 
           plot_list$Mate_fecundity + theme(legend.position="none"), 
           plot_list$MR_Female + theme(legend.position="none"), 
           plot_list$MR_Male + theme(legend.position="none"), 
           labels = "AUTO") 
# Put together plot_combined and the legend using cowplot
plot_grid(figure_2, legend_2, ncol = 1, rel_heights = c(1, .05)) 
ggsave(plot_grid(figure_2, legend_2, ncol = 1, rel_heights = c(1, .05)),
       height = 8.5,
       width = 8.5,
       file = "Output/Figure_2.jpg")

# Figure 3

# Make a combined plot without legends
figure_3 = plot_grid(plot_list$P1 + theme(legend.position="none"), 
           plot_list$P2 + theme(legend.position="none"), 
           plot_list$ML_female + theme(legend.position="none"), 
           plot_list$ML_male + theme(legend.position="none"), 
           labels = "AUTO")
# Put together plot_combined and the legend using cowplot
plot_grid(figure_3, legend_3, ncol = 1, rel_heights = c(1, .05))
ggsave(plot_grid(figure_3, legend_3, ncol = 1, rel_heights = c(1, .05)),
       height = 8.5,
       width = 8.5,
       file = "Output/Figure_3.jpg")

# Figure 4

# Make a combined plot without legends
figure_4 = plot_grid(plot_list$CD_female + theme(legend.position="none"), 
           plot_list$CD_male + theme(legend.position="none"), 
           plot_list$RML_female + theme(legend.position="none"), 
           plot_list$RML_male + theme(legend.position="none"), 
           labels = "AUTO")
# Put together plot_combined and the legend using cowplot
plot_grid(figure_4, legend_3, ncol = 1, rel_heights = c(1, .05))
ggsave(plot_grid(figure_4, legend_3, ncol = 1, rel_heights = c(1, .05)),
       height = 8.5,
       width = 8.5,
       file = "Output/Figure_4.jpg")

```
# Trade-offs between male pre- and post-copulatry traits (Table 1)

We measure additive genetic correlations between the following male traits (variable codes in brackets):

1. Sperm defence ability (p1)
2. Sperm offence ability (p2)
3. Latency to mating with virgin females (ml)
4. Latency to mating with mated females (rml)
5. Fecundity of LH females held with focal males at male biased sex ratio (matefec_mb)
6. Fecundity of LH females held with focal males at female biased sex ratio (matefec_fb)

```{r trade_offs}

# Select data for traits
d_tradeoff = d[d$Trait=="P1"|d$Trait=="P2"|d$Trait=="ML_male"|d$Trait=="RML_male"|d$Trait=="Mate_fecundity_Male_biased"|d$Trait=="Mate_fecundity_Female_biased",]

# Standardise the trait data
d_tradeoff$Measurement = scale_by(Measurement~Trait, d_tradeoff)

prior_tradeoff = list(R = list(V = diag(6), nu = 0.002), G = list(G1 = list(V = diag(6), nu = 2, alpha.mu = c(0,0,0,0,0,0), alpha.V = diag(6)*1000), G2 = list(V = 1, nu = 0.002)))
    
fit_tradeoff = MCMCglmm(Measurement ~  Trait -1, random = ~us(Trait):Family + Replicate:Trait, rcov = ~idh(Trait):units, family = "gaussian", nitt = 500000, burnin = 50000, thin=300, data = d_tradeoff, prior = prior_tradeoff, pr=TRUE, verbose = FALSE)

# Additive genetic variances for the traits

va_p1 = 2*fit_tradeoff$VCV[,"TraitP1:TraitP1.Family"]
va_p2 = 2*fit_tradeoff$VCV[,"TraitP2:TraitP2.Family"]
va_ml = 2*fit_tradeoff$VCV[,"TraitML_male:TraitML_male.Family"]    
va_rml = 2*fit_tradeoff$VCV[,"TraitRML_male:TraitRML_male.Family"]
va_matefec_mb = 2*fit_tradeoff$VCV[,"TraitMate_fecundity_Male_biased:TraitMate_fecundity_Male_biased.Family"]    
va_matefec_fb = 2*fit_tradeoff$VCV[,"TraitMate_fecundity_Female_biased:TraitMate_fecundity_Female_biased.Family"]    

# Additive genetic correlations

cor_p1_p2 = 2*fit_tradeoff$VCV[,"TraitP1:TraitP2.Family"]/(sqrt(va_p1*va_p2))
cor_p1_ml = 2*fit_tradeoff$VCV[,"TraitP1:TraitML_male.Family"]/(sqrt(va_p1*va_ml))
cor_p1_rml = 2*fit_tradeoff$VCV[,"TraitP1:TraitRML_male.Family"]/(sqrt(va_p1*va_rml))
cor_p1_matefec_mb = 2*fit_tradeoff$VCV[,"TraitP1:TraitMate_fecundity_Male_biased.Family"]/(sqrt(va_p1*va_matefec_mb))
cor_p1_matefec_fb = 2*fit_tradeoff$VCV[,"TraitP1:TraitMate_fecundity_Female_biased.Family"]/(sqrt(va_p1*va_matefec_fb))


cor_p2_ml = 2*fit_tradeoff$VCV[,"TraitP2:TraitML_male.Family"]/(sqrt(va_p2*va_ml))
cor_p2_rml = 2*fit_tradeoff$VCV[,"TraitP2:TraitRML_male.Family"]/(sqrt(va_p2*va_rml))
cor_p2_matefec_mb = 2*fit_tradeoff$VCV[,"TraitP2:TraitMate_fecundity_Male_biased.Family"]/(sqrt(va_p2*va_matefec_mb))
cor_p2_matefec_fb = 2*fit_tradeoff$VCV[,"TraitP2:TraitMate_fecundity_Female_biased.Family"]/(sqrt(va_p2*va_matefec_fb))

cor_ml_rml = 2*fit_tradeoff$VCV[,"TraitML_male:TraitRML_male.Family"]/(sqrt(va_ml*va_rml))
cor_ml_matefec_mb = 2*fit_tradeoff$VCV[,"TraitML_male:TraitMate_fecundity_Male_biased.Family"]/(sqrt(va_ml*va_matefec_mb))
cor_ml_matefec_fb = 2*fit_tradeoff$VCV[,"TraitML_male:TraitMate_fecundity_Female_biased.Family"]/(sqrt(va_ml*va_matefec_fb))

cor_rml_matefec_mb = 2*fit_tradeoff$VCV[,"TraitRML_male:TraitMate_fecundity_Male_biased.Family"]/(sqrt(va_rml*va_matefec_mb))
cor_rml_matefec_fb = 2*fit_tradeoff$VCV[,"TraitRML_male:TraitMate_fecundity_Female_biased.Family"]/(sqrt(va_rml*va_matefec_fb))

cor_matefec_mb_matefec_fb = 2*fit_tradeoff$VCV[,"TraitMate_fecundity_Male_biased:TraitMate_fecundity_Female_biased.Family"]/(sqrt(va_matefec_mb*va_matefec_fb))

if(diagnostics){
  plot(cor_p1_p2)
  print(autocorr(cor_p1_p2))
  
  plot(cor_p1_ml)
  print(autocorr(cor_p1_ml))
  
  plot(cor_p1_rml)
  print(autocorr(cor_p1_rml))
  
  plot(cor_p1_matefec_mb)
  print(autocorr(cor_p1_matefec_mb))
  
  plot(cor_p1_matefec_fb)
  print(autocorr(cor_p1_matefec_fb))
  
  plot(cor_p2_ml)
  print(autocorr(cor_p2_ml))
  
  plot(cor_p2_rml)
  print(autocorr(cor_p2_rml))
  
  plot(cor_p2_matefec_mb)
  print(autocorr(cor_p2_matefec_mb))
  
  plot(cor_p2_matefec_fb)
  print(autocorr(cor_p2_matefec_fb))
  
  plot(cor_ml_rml)
  print(autocorr(cor_ml_rml))
  
  plot(cor_ml_matefec_mb)
  autocorr(cor_ml_matefec_mb)
  
  plot(cor_ml_matefec_fb)
  print(autocorr(cor_ml_matefec_fb))
  
  plot(cor_rml_matefec_mb)
  print(autocorr(cor_rml_matefec_mb))
  
  plot(cor_rml_matefec_fb)
  print(autocorr(cor_rml_matefec_fb))
  
  plot(cor_matefec_mb_matefec_fb)
  print(autocorr(cor_matefec_mb_matefec_fb))
  
  
}


# Create Table 1

trade_off_table = matrix(NA, nrow = 6, ncol = 6)

trade_off_table[1,] = c("", "Sperm offence ability (P2)", "Latency to first mating", "Latency to mating with mated females", "Fecundity of mates (Male biased)", "Fecundity of mates (Female biased)")

trade_off_table[2, ] = c("Sperm defence ability (P1)", 
                         paste(round(mean(cor_p1_p2), 4),  " [", round(HPDinterval(cor_p1_p2)[1], 4), ",", round(HPDinterval(cor_p1_p2)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p1_ml), 4),  " [", round(HPDinterval(cor_p1_ml)[1], 4), ",", round(HPDinterval(cor_p1_ml)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p1_rml), 4),  " [", round(HPDinterval(cor_p1_rml)[1], 4), ",", round(HPDinterval(cor_p1_rml)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p1_matefec_mb), 4),  " [", round(HPDinterval(cor_p1_matefec_mb)[1], 4), ",", round(HPDinterval(cor_p1_matefec_mb)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p1_matefec_fb), 4),  " [", round(HPDinterval(cor_p1_matefec_fb)[1], 4), ",", round(HPDinterval(cor_p1_matefec_fb)[2], 4), "]", sep=""))

trade_off_table[3, ] = c("Sperm offence ability (P2)", 
                         "",
                         paste(round(mean(cor_p2_ml), 4),  " [", round(HPDinterval(cor_p2_ml)[1], 4), ",", round(HPDinterval(cor_p2_ml)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p2_rml), 4),  " [", round(HPDinterval(cor_p2_rml)[1], 4), ",", round(HPDinterval(cor_p2_rml)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p2_matefec_mb), 4),  " [", round(HPDinterval(cor_p2_matefec_mb)[1], 4), ",", round(HPDinterval(cor_p2_matefec_mb)[2], 4), "]", sep=""),
                         paste(round(mean(cor_p2_matefec_fb), 4),  " [", round(HPDinterval(cor_p2_matefec_fb)[1], 4), ",", round(HPDinterval(cor_p2_matefec_fb)[2], 4), "]", sep=""))

trade_off_table[4, ] = c("Latency to first mating", 
                         "",
                         "",
                         paste(round(mean(cor_ml_rml), 4),  " [", round(HPDinterval(cor_ml_rml)[1], 4), ",", round(HPDinterval(cor_ml_rml)[2], 4), "]", sep=""),
                         paste(round(mean(cor_ml_matefec_mb), 4),  " [", round(HPDinterval(cor_ml_matefec_mb)[1], 4), ",", round(HPDinterval(cor_ml_matefec_mb)[2], 4), "]", sep=""),
                         paste(round(mean(cor_ml_matefec_fb), 4),  " [", round(HPDinterval(cor_ml_matefec_fb)[1], 4), ",", round(HPDinterval(cor_ml_matefec_fb)[2], 4), "]", sep=""))


trade_off_table[5, ] = c("Latency to mating with mated females", 
                         "",
                         "",
                         "",
                         paste(round(mean(cor_rml_matefec_mb), 4),  " [", round(HPDinterval(cor_rml_matefec_mb)[1], 4), ",", round(HPDinterval(cor_rml_matefec_mb)[2], 4), "]", sep=""),
                         paste(round(mean(cor_rml_matefec_fb), 4),  " [", round(HPDinterval(cor_rml_matefec_fb)[1], 4), ",", round(HPDinterval(cor_rml_matefec_fb)[2], 4), "]", sep=""))


trade_off_table[6, ] = c("Fecundity of mates (Male biased)", 
                         "",
                         "",
                         "",
                         "",
                         paste(round(mean(cor_matefec_mb_matefec_fb), 4),  " [", round(HPDinterval(cor_matefec_mb_matefec_fb)[1], 4), ",", round(HPDinterval(cor_matefec_mb_matefec_fb)[2], 4), "]", sep=""))

# save file

write.csv(noquote(trade_off_table), file = "Output/trade_off_output.csv", row.names = FALSE, col.names = FALSE)
print(data.frame(trade_off_table))


```

# Analysing female fecundity (Table S5)

Here we investigate the effect of sex ratio on the fecundity of females from the LH population.

```{r fecundity}
d_fec = d[d$Trait=="Mate_fecundity_Male_biased"|d$Trait=="Mate_fecundity_Female_biased",]
d_fec$Sex.Ratio = ifelse(d_fec$Trait=="Mate_fecundity_Male_biased", "Male biased", "Female biased")

fit_fec = lmer(Measurement ~ Sex.Ratio + (1|Family:Sex.Ratio) + (1|Replicate), d_fec)
summary(fit_fec)

#################
### Figure S3 ###
#################

p_fec = ggplot(d_fec, aes(y = Measurement, x = Sex.Ratio))
p_fec = p_fec + theme_bw() + geom_boxplot() + theme(text = element_text(size = 15)) + labs(x = "Sex ratio", y = "Fecundity of LH females")
p_fec
ggsave(p_fec, 
       height = 6.5,
       width = 6.5,
       file = "Output/Figure_S3.jpg")
```

